  set more off
  mat drop _all
  set seed 19112014

*===================================================================*
*   BIHAR EVALUATION OF SOCIAL FRANCHISING AND TELEMEDICINE (BEST)
*                           DD ANALYSIS
*====================================================================*

*********************************************************************
* Table 2: DD Analysis
* Format Mean(SD) for each group, Difference (95% CI) and p-value
*********************************************************************

*====================================================================*
* Diarrhea
*====================================================================*

* Append Datasets
  use "$prodata1\providers_interview_1st", clear
  merge 1:1 prov_id wave using "$prodata1\vignette_diarrhea", keepus(case sample4 nhistory dx_type tx_any tx_detail tx_type vig_theta_mle_q12_std vig_item? vig_item??)
  drop _merge 

  append using "$prodata2\providers_interview_2nd"
  merge 1:1 prov_id wave using "$prodata2\vignette_diarrhea", update keepus(case sample4 nhistory dx_type tx_any tx_detail tx_type vig_theta_mle_q12_std vig_theta_mle_q12_stdb vig_item? vig_item??)
  drop _merge 

* Checking district in wave 2
  assert cluster == 150 if district==2 & wave==2
  assert cluster == 322 if district==5 & wave==2

* Sample (see do file 3_samples)
  sort wave prov_id
  merge 1:1 wave prov_id using "$prodata\sample_by_wave"
  keep if _merge==3
  drop _merge 
  assert _N==810

* Sample selection 2
/*NOTE: we need to exclude WHP providers*/
  gen     regsample = 1
  replace regsample = 0 if WHP_oversampled==1

* Post variable
  gen post = (wave==2)

* Interaction 
  gen treat_post = treat*post

* Percentage of selected items asked
  egen pq_vig_d=rmean(vig_item1-vig_item10 vig_item12-vig_item14)

* Unnecessary and harmful treatments 
  gen tx_harm=(tx_detail>=3 & tx_detail<=5 | tx_detail==7)

* ORS
  gen     ORS = (tx_detail<=5)
  lab var ORS "=1 if ORS given"

* Knowledge Score
  gen     knowscore = vig_theta_mle_q12_std if wave==1
  replace knowscore = vig_theta_mle_q12_stdb if wave==2

* Dependent variable
  local dep = "pq_vig_d dx_type ORS tx_type tx_harm"

* Regressors
  local policy "treat post treat_post"
  local reg1 "`policy'"
  local reg2 "`policy' age experience mqualif knowscore"
  local reg3 "`policy' age experience mqualif knowscore work_hs npatient camps fac_gov niindex"

* Matrix to store results
  mat table2 = J(11,11,.)
  mat colnames table2 = N N_pre mean_pre N_post mean_post beta se ci_l ci_u pvalue adj_pvalue


* Regressions
  local i 1
  foreach var of local dep {

    * Mean value of dependent variable 
      sum `var' if regsample==1 & post==0
      local npre = r(sum)
      local mean_pre = string(round(r(mean)*100, .1))

      sum `var' if regsample==1 & post==1
      local npost = r(sum)
      local mean_post = string(round(r(mean)*100, .1))


      forvalues r=1/3 {

      * Regressions
        areg `var' `reg`r'' if regsample==1, absorb(district) cluster(cluster)

      * Fill matrix
        if `r'==3 {
          mat table2[`i',1] = e(N)
          mat table2[`i',2] = `npre'
          mat table2[`i',3] = `mean_pre'
          mat table2[`i',4] = `npost'
          mat table2[`i',5] = `mean_post'
          mat table2[`i',6] = _b[treat_post]
          mat table2[`i',7] = _se[treat_post]
          local low_ci  = _b[treat_post] - invttail(e(df_r),0.025)*_se[treat_post]
          local high_ci = _b[treat_post] + invttail(e(df_r),0.025)*_se[treat_post]
          mat table2[`i',8] = `low_ci'
          mat table2[`i',9] = `high_ci'
          local t = _b[treat_post]/_se[treat_post]
          mat table2[`i',10] = 2*ttail(e(df_r),abs(`t'))
        }

    * Save results
      if "`dep'"=="pq_vig_d" & `r'==1 local app = "replace"
      else local app = "append"
      outreg2 using "$results\table2_diarrhea_vignette.xls", keep(treat_post) ctitle(`var') bdec(3) sdec(3) ///
      nocons stats(coef ci_low ci_high pval) noparen addnote("District FE") `app'

    }
  local i = `i'+1
  }



* Check: cluster fixed effect instead of treat
  local reg "age experience mqualif knowscore work_hs npatient camps fac_gov niindex"
  mat table2_clusterFE = J(11,11,.)
  mat colnames table2_clusterFE = N N_pre mean_pre N_post mean_post beta se ci_l ci_u pvalue adj_pvalue
  local i 1
  foreach var of local dep {

    * Mean value of dependent variable 
      sum `var' if regsample==1 & post==0
      local npre = r(sum)
      local mean_pre = string(round(r(mean)*100, .1))

      sum `var' if regsample==1 & post==1
      local npost = r(sum)
      local mean_post = string(round(r(mean)*100, .1))

    * Regression
      areg `var' `reg' i.cluster post treat_post if regsample==1, absorb(district) cluster(cluster)


    * Save results
      outreg2 using "$results\table2_diarrhea_vignette.xls", keep(treat_post) ctitle("cluster FE",`var') bdec(3) sdec(3) ///
      stats(coef ci_low ci_high pval) noparen nocons append


    * Fill matrix
      mat table2_clusterFE[`i',1] = e(N)
      mat table2_clusterFE[`i',2] = `npre'
      mat table2_clusterFE[`i',3] = `mean_pre'
      mat table2_clusterFE[`i',4] = `npost'
      mat table2_clusterFE[`i',5] = `mean_post'
      mat table2_clusterFE[`i',6] = _b[treat_post]
      mat table2_clusterFE[`i',7] = _se[treat_post]
      local low_ci  = _b[treat_post] - invttail(e(df_r),0.025)*_se[treat_post]
      local high_ci = _b[treat_post] + invttail(e(df_r),0.025)*_se[treat_post]
      mat table2_clusterFE[`i',8] = `low_ci'
      mat table2_clusterFE[`i',9] = `high_ci'
      local t = _b[treat_post]/_se[treat_post]
      mat table2_clusterFE[`i',10] = 2*ttail(e(df_r),abs(`t'))

      local i = `i'+1
  }


* Erase txt files created by outreg
  erase "$results\table2_diarrhea_vignette.txt"



*====================================================================*
* Pneumonia
*====================================================================*

* Append Datasets
  use "$prodata1\providers_interview_1st", clear
  merge 1:1 prov_id wave using "$prodata1\vignette_pneumonia", keepus(case sample4 nhistory dx_type tx_any tx_detail tx_type vig_theta_mle_q12_std vig_item? vig_item??)
  drop _merge

  append using "$prodata2\providers_interview_2nd"
  merge 1:1 prov_id wave using "$prodata2\vignette_pneumonia", update keepus(case sample4 nhistory dx_type tx_any tx_detail tx_type vig_theta_mle_q12_std vig_theta_mle_q12_stdb vig_item? vig_item??)
  drop _merge

* Sample
  sort wave prov_id
  merge 1:1 wave prov_id using "$prodata\sample_by_wave"
  keep if _merge==3
  drop _merge 
  assert _N==810

* Sample selection 2
/*NOTE: we need to exclude WHP providers*/
  gen     regsample = 1
  replace regsample = 0 if WHP_oversampled==1

* Post variable
  gen post = (wave==2)

* Interaction 
  gen treat_post = treat*post

* Percentage of selected items asked
  egen pq_vig_p=rmean(vig_item1-vig_item9)

* Unnecessary and harmful treatments 
  gen tx_harm=(tx_detail>=3 & tx_detail<=5)

* Harmful or non treatment
  gen tx_harm_no=(tx_detail>=3 & tx_detail<=6)

* Antibiotics
  gen antibiotics=(tx_detail<=4)

* Knowledge Score
/*NOTE: For the second round we need the score without considereing the oversampled WHP*/
  gen     knowscore = vig_theta_mle_q12_std if wave==1
  replace knowscore = vig_theta_mle_q12_stdb if wave==2

* Dependent variable
  local dep = "pq_vig_p dx_type antibiotics tx_type tx_harm tx_harm_no"
  local dep = "pq_vig_p dx_type antibiotics tx_type tx_harm"

* Regressors
  local policy "treat post treat_post"
  local reg1 "`policy'"
  local reg2 "`policy' age experience mqualif knowscore"
  local reg3 "`policy' age experience mqualif knowscore work_hs npatient camps fac_gov niindex"

* Regressions
  local i 7
  foreach var of local dep {
    forvalues r=1/3 {

    * Mean value of dependent variable 
      sum `var' if regsample==1 & post==0
      local npre = r(sum)
      local mean_pre = string(round(r(mean)*100, .1))

      sum `var' if regsample==1 & post==1
      local npost = r(sum)
      local mean_post = string(round(r(mean)*100, .1))

    * Regressions
      areg `var' `reg`r'' if regsample==1, absorb(district) cluster(cluster)

    * Fill matrix
      if `r'==3 {
        mat table2[`i',1] = e(N)
        mat table2[`i',2] = `npre'
        mat table2[`i',3] = `mean_pre'
        mat table2[`i',4] = `npost'
        mat table2[`i',5] = `mean_post'
        mat table2[`i',6] = _b[treat_post]
        mat table2[`i',7] = _se[treat_post]
        local low_ci  = _b[treat_post] - invttail(e(df_r),0.025)*_se[treat_post]
        local high_ci = _b[treat_post] + invttail(e(df_r),0.025)*_se[treat_post]
        mat table2[`i',8] = `low_ci'
        mat table2[`i',9] = `high_ci'
        local t = _b[treat_post]/_se[treat_post]
        mat table2[`i',10] = 2*ttail(e(df_r),abs(`t'))
      }


    * Save results
      if "`dep'"=="pq_vig_p" & `r'==1 local app = "replace"
      else local app = "append"
      outreg2 using "$results\table2_pneumonia_vignette.xls", keep(treat_post) ctitle(`var') bdec(3) sdec(3) ///
      nocons stats(coef ci_low ci_high pval) noparen addnote("District FE") `app'

    }
    local i = `i'+1
  }

* Check: cluster fixed effect instead of treat
  local reg "age experience mqualif knowscore work_hs npatient camps fac_gov niindex"
  local i 7
  foreach var of local dep {

    * Mean value of dependent variable 
      sum `var' if regsample==1 & post==0
      local N_pre = r(sum)
      local mean_pre = string(round(r(mean)*100, .1))

      sum `var' if regsample==1 & post==1
      local N_post = r(sum)
      local mean_post = string(round(r(mean)*100, .1))

    * Regressions
      areg `var' `reg' i.cluster post treat_post if regsample==1, absorb(district) cluster(cluster)

    * Save results
      outreg2 using "$results\table2_pneumonia_vignette.xls", keep(treat_post) ctitle("cluster FE",`var') bdec(3) sdec(3) ///
      nocons stats(coef ci_low ci_high pval) noparen append

    * Fill matrix
      mat table2_clusterFE[`i',1] = e(N)
      mat table2_clusterFE[`i',2] = `npre'
      mat table2_clusterFE[`i',3] = `mean_pre'
      mat table2_clusterFE[`i',4] = `npost'
      mat table2_clusterFE[`i',5] = `mean_post'
      mat table2_clusterFE[`i',6] = _b[treat_post]
      mat table2_clusterFE[`i',7] = _se[treat_post]
      local low_ci  = _b[treat_post] - invttail(e(df_r),0.025)*_se[treat_post]
      local high_ci = _b[treat_post] + invttail(e(df_r),0.025)*_se[treat_post]
      mat table2_clusterFE[`i',8] = `low_ci'
      mat table2_clusterFE[`i',9] = `high_ci'
      local t = _b[treat_post]/_se[treat_post]
      mat table2_clusterFE[`i',10] = 2*ttail(e(df_r),abs(`t'))

      local i = `i'+1
  }

* Erase txt files created by outreg
  erase "$results\table2_pneumonia_vignette.txt"

* Export tables
  local matrix "table2 table2_clusterFE"
  foreach mat of local matrix {
  preserve
  drop _all
  svmat2 `mat', names(col)
  foreach var of varlist beta-adj_pvalue {
    replace `var' = round(`var',.001)
  }
  outsheet using "$results\\`mat'.xls", replace
  restore
  }
